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1.  Summary 


The  Radiative  Transfer  (RT)  approach  is  widely  used  in  applications  involving  scattering 
from  layered  random  media  with  rough  interfaces  (F.T.  Ulaby,  R.K.  Moore  and  AK  Fung, 
Microwave  Remote  Sensing,  Vol.  3,  Artech  House,  1986).  Although  this  approach  involves 
approximations,  in  most  applications  they  are  not  explicitly  stated  or  well  understood.  The 
RT  approach  for  random  media  with  non-scattering  boundaries  has  been  well  studied  (M.I. 
Mishchenko,  Appl.  Optics,  41,  7114-34,  2002).  In  contrast  our  problem  has  scattering 
boundaries  which  are  randomly  rough.  In  order  to  better  understand  the  RT  approach  to 
our  problem  we  adopt  a  statistical  wave  approach  and  then  transition  to  the  RT  equations. 
The  geometry  of  our  problem  consists  of  a  multi-layer  discrete  random  medium  with  rough 
boundaries  which  are  planar  on  the  average.  The  random  medium  in  each  layer  consists  of 
a  homogeneous  background  medium  in  which  discrete  scatterers  are  randomly  distributed. 
The  statistical  characteristics  of  the  random  medium  in  each  layer  are  independent  of  each 
other  and  independent  of  the  statistics  describing  the  rough  interfaces.  The  regions  above 
and  below  the  random  medium  stack  are  homogeneous.  Using  the  Greens  functions  of  the 
problem  without  the  volumetric  fluctuations  we  represent  our  problem  as  a  system  of 
integral  equations.  Employing  the  T-matrix  description  we  first  average  with  respect  to 
volumetric  fluctuations  and  then  use  the  Twersky  approximation  to  obtain  a  system  of 
integral  equations.  We  next  average  with  respect  to  surface  fluctuations,  apply  the  weak 
surface  correlation  approximation  and  arrive  at  a  closed  system  of  integral  equations  for 
the  first  and  second  moments  of  the  electric  fields.  We  use  the  Wigner  transforms  to 
translate  the  coherence  functions  to  radiant  intensities,  which  are  the  fundamental 
quantities  in  the  RT  approach.  On  applying  the  quasi-static  field  approximation  we  hence 
arrive  at  a  system  of  equations  identical  to  those  used  in  the  RT  approach.  From  this  study 
we  leam  that  there  are  more  conditions  involved  in  the  RT  approach  than  widely  believed 
to  be  sufficient.  This  means  that  the  applicability  of  the  RT  approach  is  much  smaller  than 
what  is  generally  conceived  by  the  remote  sensing  community. 
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2.  Introduction 


The  radiative  transfer  (RT)  theory  is  widely  used  in  remote  sensing  problems 
[16,28,17,15,1,38].  Often  the  model  of  layered  random  medium  with  rough  interfaces  is 
used.  Multiple  scattering  processes  in  this  structure  are  well  represented  by  the  RT 
equations.  Although  quite  successful  in  numerous  applications  in  various  disciplines,  it  is 
known  that  the  RT  approach  involves  approximations.  Often  people  in  the  remote  sensing 
community  are  not  quite  familiar  with  the  approximations  involved  in  the  RT  approach  and 
hence  there  has  been  inappropriate  use  of  the  RT  approach  in  the  literature.  Since  the 
phenomenological  RT  theory  [7,35]  was  first  developed  for  light  scattering  in  planetary 
atmospheres  the  RT  conditions  prevalent  in  the  atmospheric  context  has  been  popularly 
identified  as  sufficient  conditions  for  employing  the  RT  theory.  However,  we  notice  that 
the  RT  theory  has  been  used  for  a  variety  of  different  problems  [25,32,26,34]  in  various 
applications  with  complex  geometries.  It  is  not  clear  whether  the  classical  conditions 
associated  with  the  RT  theory  are  sufficient  in  all  situations.  In  this  report  we  will  review 
the  approximations  involved  and  clarify  misconceptions.  In  order  to  better  understand  the 
RT  approach  we  employ  the  more  rigorous  statistical  wave  theory  to  the  problem  and 
hence  make  the  transition  to  the  RT  equations.  In  this  process  we  clarify  and  explain  the 
assumptions  or  approximations  involved  in  the  RT  approach.  By  following  this  procedure 
we  found  that  there  are  more  conditions  embedded  in  the  RT  approach  than  widely 
believed  to  be  sufficient.  For  our  study  we  have  considered  a  multi-layer  random  medium 
composed  of  discrete  scatterers  (see  Figure  1).  By  considering  several  special  cases  of  this 
general  problem  we  show  that  the  number  of  conditions  implied  in  the  RT  approach 
reduces  with  simpler  geometries.  Our  conclusions  are  not  just  for  the  case  discrete  random 
media.  Similar  conditions  apply  for  the  random  continuum  as  well.  Similar  conditions 
apply  for  the  random  continuum  as  well  [27].  There,  the  volumetric  fluctuations  were 
modeled  as  a  random  continuum.  The  random  continuum  concept  is  rather  abstract  and 
difficult  to  model  several  problems  encountered  in  applications.  Perhaps  the  most  serious 
drawback  of  the  continuum  model  is  that  it  often  impossible  to  relate  the  key  statistical 
parameters  to  observable  quantities.  The  current  report  is  devoted  to  the  remote  sensing 
applications.  A  discrete  scatterer  model  is  most  appropriate  and  realistic  for  modeling  real 
world  problems  [24].  The  analysis  used  in  the  discrete  random  media  is  quite  different  from 
the  continuum  case.  Furthermore,  the  conditions  and  assumptions  embedded  in  the  RT 
approach  are  also  different  in  the  discrete  model  case. 


The  report  is  organized  as  follows.  First  we  describe  the  geometry  of  the  problem.  Next  we 
give  the  radiative  transfer  approach  to  the  problem.  The  next  section  is  on  the  statistical 
wave  approach  to  the  problem.  This  occupies  the  major  part  of  the  report.  It  deals  with  the 
derivation  and  analysis  of  the  first  and  second  moments  of  the  wave  functions.  A  transition 
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is  next  made  to  RT  equations.  Next  we  enumerate  and  discuss  the  conditions  implied  in  the 
RT  approach. 
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3.  Methods,  Assumptions,  and  Procedures 

3.1  Description  of  the  Problem 

The  geometry  of  the  problem  is  shown  in  Figure  1.  We  have  an  /V-laycr  random  medium 
stack  with  rough  interfaces  which,  on  the  average,  are  parallel  planes.  Let  e  be  the 

permittivity  of  the  background  medium,  and  let  s  be  the  permittivity  of  the  scatterers  in 

the y'-th  layer.  The  location  and  orientation  of  the  scatterers  are  random  functions 
characterizing  the  fluctuations.  On  the  average  the  problem  is  translationally  invariant  and 
isotropic  in  the  azimuth.  We  assume  that  the  volumetric  fluctuations  in  each  layer  are 
statistically  independent  of  each  other.  Let  N}  be  the  number  of  scatterers,  and  let  p.  be  the 

density  of  the  scatterers  (number  of  scatterers  per  unit  volume)  of  the  /-the  layer.  The 
permeability  of  all  the  layers  is  that  of  free  space.  The  randomly  rough  interfaces  are  given 
as  z  =  z  .  +  £,(r  ) .  Then  .  arc  zero-mean  isotropic  stationary  random  processes 

independent  of  volumetric  fluctuations  of  the  problem.  Let  zo  =  0,  and  let  d}  be  the 
thickness  of  the  y'-th  layer.  Let  zo  =  0,  and  let  d}  be  the  thickness  of  the  y'-th  layer.  The  media 
above  and  below  the  stack  are  homogeneous  with  parameters  <w ,  ko,  and  sN+l ,  k\+\, 

respectively.  This  system  is  excited  by  a  monochromatic  electromagnetic  plane  wave  and 
we  are  interested  in  formulating  the  resulting  multiple  scattering  process. 
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Figure  1.  Geometry  of  the  Problem 


3.2  Radiative  Transfer  Approach 

Multiple  scattering  in  a  complex  environment  is  well  described  by  the  radiative  transfer 
theory.  This  theory  is  not  only  conceptually  simple  but  also  very  efficient.  The  fundamental 
quantity  here  is  the  specific  intensity  I  which  is  governed  by  the  following  equation 
[7,35,11] 

s-VI(r,5)  +  yI(r,i)  =  |P(i,5')I(r,i')(if2'  (0 

One  may  regard  this  equation  as  a  statement  of  conservation  of  energy  density  I  which  is  a 
phase-space  quantity  at  position  r  and  direction  s  .  y  is  the  extinction  matrix  which  is  a 

measure  of  loss  of  energy  due  to  scattering  in  other  directions.  P  is  the  phase  matrix 
representing  increase  in  energy  density  due  to  scattering  from  neighbouring  elements.  Q  is 
the  solid  angle  subtended  by  s  .  Given  the  statistical  characteristics  of  the  medium  one  can 
calculate  the  phase  function  using  the  single  scattering  theory  for  elements  that  constitute 
the  random  medium  of  the  layer  [39,12,6,23].  The  extinction  matrix  is  hence  calculated 
using  the  optical  theorem.  The  specific  intensity  in  each  layer  is  governed  by  an  equation 
similar  to  (1).  Since  our  layer  problem  has  translational  invariance  in  azimuth  the  RT 
equation  for  the  m- th  layer  takes  the  following  form, 

cos0Tz l-(Zii)+^i-(z,S)  =  J0.  <2) 


where  the  subscript  m  denotes  that  the  quantity  corresponds  to  that  of  the  m-th  layer  and  0 
is  the  elevation  angle  of  s.  This  set  of  RT  equations  is  complemented  by  a  set  of  boundary 
conditions  which  is  in  turn  based  on  energy  conservation  considerations.  To  be  more 
precise,  we  impose  the  condition  that  the  energy  flux  density  at  each  interface  is  conserved. 
This  leads  to  the  following  boundary  conditions  on  the  m-th  interface 


i;  (z.  J) = 1  (s.s’)  i;  + 1  x.,„.  a,  s')  i 


(3) 


The  boundary  conditions  on  the  (m  - 1)  -th  interface  are  given  as 


=  JX.u.ftO  C(z„-,  s' W  +  1^,(7..,  i1 yci' 


(4) 
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where  Rmn  and  Xmn  are  the  local  reflection  and  transmission  Mueller  matrices.  To  be  more 

specific,  R(ot  represents  the  reflection  matrix  of  waves  incident  from  medium  n  on  the 

interface  separating  medium  m  and  medium  n.  The  superscripts  it  and  d  indicate  whether 
the  intensity  corresponds  to  a  wave  travelling  upwards  or  downwards.  These  Mueller 
matrices  are  often  calculated  using  some  asymptotic  theory  such  as  the  Kirchhoff 
approximation  [41,5,38,37].  The  integration  in  these  expressions  are  over  a  solid  angle 
(hemisphere)  corresponding  to  s ' .  Suppose  we  have  a  time-harmonic  electromagnetic  plane 
wave  incident  on  this  stack  from  above.  Then  the  downward  travelling  intensity  in  Region 
0  is 


I d0(z,s)  =  Boc>(cos0o  -  cos  -$)  (5) 

where  B0  is  the  intensity  of  the  incident  plane  wave  and  { 8j ,  (pt }  describes  its  direction. 
Since  there  is  no  source  or  scatterer  in  Region  N+ 1, 

l^+1(z,s)  =  0 

Notice  again  that  these  boundary  conditions  represent  conservation  of  intensity  at  the 
interfaces.  We  should  point  out  that  the  radiative  transfer  approach  as  applied  to  a 
particular  problem  is  only  a  model  based  on  certain  assumptions.  Since  the  RT  theory  is 
used  in  a  variety  of  applications,  the  particular  assumptions  involved  are  described  in  terms 
of  different  tenninologies,  specific  to  the  discipline  where  it  is  used.  One  good  way  to 
understand  in  more  general  terms  the  RT  approach  and  the  underlying  assumptions  is  to 
compare  it  with  the  more  rigorous  wave  approach.  For  the  case  of  an  unbounded  random 
medium  this  kind  of  study  was  carried  out  in  the  1970s  [3,2],  From  that  study  we  learn  that 
the  radiative  transfer  theory  can  be  applied  under  the  following  conditions: 

1 .  Quasi-stationary  field  approximation 

2.  Sparse  distribution 

3.  Statistical  homogeneity  of  the  medium  fluctuations 

These  are  the  well-known  conditions  that  we  associate  with  the  RT  approach.  However,  our 
problem  has  bounded  structures  and,  further,  they  are  randomly  rough.  The  question  is  this: 
are  the  above  conditions  sufficient  to  apply  the  RT  approach  for  our  problem?  This  is  the 
motivation  for  this  report.  We  follow  the  wave  approach  to  this  problem,  derive  the 
equations  for  the  intensities,  and  hence  make  the  transition  to  the  RT  equations.  This 
procedure  enables  us  to  better  understand  the  necessary  conditions  for  using  the  RT 
approach  for  our  problem. 
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3.3  Statistical  Wave  Approach 

The  following  are  the  equations  that  govern  the  waves  in  the  layer  structure: 


V  x  V  x  E  j  -  k2j E j  =  v;.E  j 

j=h-,N 

(6) 

where 

v;=ZQji(r) 

1=1 

(7a) 

r  €  Vji 


otherwise 


(7b) 


k)s  =  (O ljU£js  k)  =  (Q1  /U£j  (?c) 

V -  is  the  domain  of  the  z'-th  scatterer  and  N  is  the  number  of  scatterers  in  layer  /.  Thus  v 

Jl  J  J  J  J 

represents  the  volumetric  fluctuations  in  Region  j.  For  the  homogeneous  regions  above  and 


below  we  have 

VxVxE0-k02E0  =  0  (8a) 

VxVxE*+1-*’+1E*+1=0  (8b) 

The  boundary  conditions  at  the y'-th  interface  are 

n  x  Ej  (r±  ,£j)  =  fix  EJ+l  (r± ,  Cj)  (9a) 

hxVxE.(r±,4'y.)  =  nxVxEi+1(r1,4'.)  (9b) 


where  h  is  the  unit  vector  normal  to  the  y-th  interface  with  nonnal  pointing  into  the 
medium  j.  This  system  is  complemented  by  the  radiation  conditions  well  away  from  the 
stack.  We  assume  that  we  know  the  solution  to  the  problem  without  volumetric 
fluctuations,  and  denote  it  as  E  .  Let  the  Green's  functions  to  this  problem  be  denoted  as 

G.. .  These  Green's  functions  are  governed  by  the  following  set  of  equations: 

V  x  V  x  e Jt (r,  r')  - k]& (r,  r’)  =  %<5(r  -  r')  < 1 
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(10b) 


n  x  GJk  (r± ,  Cj  ;r')  =  nx  Gy+1)t  (r± ,  £ ) ;  r') 
n  x  V  x  tjk  (r± ,  £j ;  r')  =  n  x  V  x  6,  /+1),  (r± ,  £j ;  r') 


(10c) 


where  I  is  unit  dyad.  The  boundary  conditions  above  are  for  the  j- th  interface.  Similarly 
on  the  (/- 1  )-th  interface  we  have  the  following  boundary  conditions. 


n x Gjk (r  ,£J_l;r')  =  nxG( M), (rx , £M ; r ') 


{j-\)k  V  _L5  j- P 


(11a) 


nxVxG,  (r± ,  C_,;r')  =  nxVxG.  M)t  (rx ,  ^  ;  r') 


(lib) 


Using  these  Green's  functions  and  the  radiation  conditions  the  wave  functions  can  be 
represented  as 


N 

E,(r)  =  E,(r)  +  X( 

k= 1  ‘ 


Q, 


G.,(r,rX(r')E,(r')Jr' 


j  =  0,1,...,  TV +  1 


(12) 


where  Qk  -  {r [;£k  <z  <  C/:  , } .  Note  that  v0  =  vv+]  =  0 .  In  order  to  carry  out  multiple 
scattering  analysis  with  a  distribution  of  discrete  scatterers  it  is  convenient  to  employ  the 
concept  of  transition  operator  T  [15,42,40].  Suppose  we  know  the  electric  field  E;  incident 
on  the  /-th  scatterer.  We  introduce  the  transition  operator  T/  such  that  the  electric  field 
scattered  by  the  /-th  scatterer  is  given  as  TX  .  Using  this  concept  (12)  may  be  expressed  as 


N  Nk 


E.(r)  =  E.(r)  +  ZEL  dr'\n  dr"Gjk  (r,  r')Tj  (r',  r")E'  (rff)  j  =  0, 1, . . . ,  N  + 1  (]3) 


k= 1  1=1 


Note  that  Tlk  depends  only  on  the  /-th  scatterer.  To  proceed  further  with  the  multiple 
scattering  analysis  it  is  expedient  to  use  symbolic  representation  of  (13). 

E,=k+e,<-Tx  (14) 


First  we  average  (14)  w.r.t.  volumetric  fluctuations  to  get 
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where  the  subscript  v  denotes  volumetric  averaging.  Since  there  are  Nk  scatterers  in  the  k- 
th  layer 


^  —  J  rYtEJkp(jcl ,  r2 , . . . ,  ,  Sj ,  s2 , . . . ,  )  drds 


(16) 


where  p  is  the  j oint  probability  density  function  of  finding  the  scatterers  at  r, ,  r2 , . . . ,  rv 
with  orientations  Sj ,  s2 , . . . ,  s N  .  We  assume  that  the  positions  and  orientations  are 
independent  of  each  other.  In  other  words 


p(?\ 1 1*2’ •  •  •  ’ rv„ )  piy j , r2,. . . , rN  )/?(Sj, s2, . . . , s„  ) 


(17) 


Furthermore,  assume  that  the  orientation  of  the  particle  at  position  r(  is  independent  of  the 
orientation  of  all  other  particles,  which  means 


P(  S1,S2,' 


s)=Z^(si) 


We  next  express  the  joint  position  probability  density  function  as 
p(rl,r2,...,rNk)  =  p(rvr2,...,r;,...,rNk\r,\)p(ri) 


(18) 


(19) 


where  p(rt ,  r2 , . . . ,  r/ , . . . ,  rv.  |  r, )  is  the  conditional  pdf  of  finding  scatterers  at  r, ,  r2 , . . . ,  rv 

by  fixing  the  /-th  scatterer  is  at  r, .  The  prime  in  r/  denotes  that  r,  should  be  omitted  in  the 

argument  list  of  the  conditional  probability  density  function.  Substituting  this  relation  in 
(16)  we  obtain 

(TX)  =  (T,'{E;)j  <20) 

where  \  denotes  conditional  average  with  scatterer  I  fixed  at  r, .  If  W  is  large  and  the 
distance  between  scatterers  is  large  then  we  can  approximate 

(e;)'d(e;.)  <21> 

This  is  called  the  Foldy’s  approximation  and  is  applicable  for  sparse  media.  Under  this 
approximation  (15)  becomes 
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(E,)>DE,+AGJl{Tt)(E1> 


(22) 


where  pk  is  the  density  of  the  scatterers  in  layer  k.  Now  operate  the  above  equation  by 
V  x  V  x  I  -  kj  I  to  obtain 


vxVx(E,y-*;(E,y=pyT,)(E 


(23) 


Next  average  this  over  surface  fluctuations  to  get 


VxVx 


(E/>, 


E, )  =P 


f)(E, 


(24) 


Note  that  the  transition  operators  are  independent  of  surface  fluctuations.  From  this  we  see 
that 


Vx Vx(E .)  -kUE,)  =0  j  =  0,N  +  l  (25) 

which  means  that  the  coherent  propagation  constants  in  the  regions  above  and  below  the 
layered  stack  are  unaffected  by  the  fluctuations  of  the  problem.  However,  they  indeed  get 
modified  in  the  layered  stack  region.  On  writing  (24)  as 


{Vx  Vxl  -k*  (t,)}(e,.)^  =  0 


(26) 


we  infer  that  x}  =  -yjk/  +  P,  (T)  }  represents  the  mean  propagation  constant,  in  operator 
form,  for  coherent  waves  in  layer  j. 


Since  the  problem  is  statistically  homogeneous  in  azimuth  the  mean  fields  in  our  system 
have  the  following  form: 


(Et(r))v,  =exP0'k±rr) 

{ (kjj  )p+  exp  \iqpjz\  +  fif  (k±i  )p^  exp  [-igjz]}  j  =  1, 2, . . . ,  N 


(27) 


{K(r))vs  =  exp(/kii  -r)|p0  exp[-^.z]  +  i?p(k±!.)p;  exp [^,-z]} 
(Ev+i(r))vs  =  exp(zk_u  •r)Xp(k±,.)p^+1  exp[-/k/;+l)z!z] 
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(28) 

(29) 


where  the  superscript  p  stands  for  the  polarization,  either  horizontal  or  vertical,  p  is  the 
unit  vector  representing  polarization,  qj  is  the  z-component  of  .  The  subscript  i  is  used  to 

indicate  that  the  wave  vector  is  in  the  incident  direction.  R  and  X  denote  respectively  the 
mean  reflection  and  transmission  coefficient  of  the  stack.  A  .  and  B  denote  respectively  the 

mean  coefficients  of  up-going  and  down-going  waves  in  the y'-th  layer. 


Based  on  this  we  can  formulate  the  waves  averaged  w.r.t.  volumetric  fluctuations  as 


(Ew), 


=^jlAiexpo'ki'r){'4r(ki.kil)9/ e',";+5r(ki>ki/>«,7 e 

=  exp(/k±  •r)exp[-/k0ziz]p“  +~rj  exp(/k±  •r)/?w(k1,k1()y(;  exp[ik0zz]dk 


(30) 

(31) 


and 

(®v+i(r))v  =  exp0'k± -r)Xpq(k±,k±i)qN+1  exp[-ik(NU)zz]dk±  (32) 


where  Af.  ,BrR  and  X  are  now  integral  operators  representing  scattering  from  rough 

interfaces.  The  boundary  conditions  associated  with  the  above  equations  at  the  j- th 
interface  are 

nx(E.(r±,^))v=nx(E.+1(ri,^))v 
and 

nxVx(E.  (r± ,  Cj  =  n  x  V  x  (E/+l  (r± ,  Cj))v 

The  above  system  may  be  solved  either  numerically  or  by  any  one  of  the  asymptotic 
methods  available  in  rough  surface  scattering  theory  [5,4,41]  to  evaluate  the  mean 
coefficients  that  appear  in  (27)-(29). 

We  proceed  now  to  the  analysis  of  the  second  moments,  by  starting  with  (12).  For 
convenience  we  write  it  in  symbolic  form  as 


7=1,2, 


,N 


(33) 


7  =1,2,..., 


A  \ 
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(35) 


N 


E/=%+ZG,.-TX 


k=l 


We  take  the  tensor  product  of  this  equation  with  its  complex  conjugate  and  average  w.r.t. 
volumetric  fluctuations  and  obtain 

N  N  Nk  Nt. 

k= 1  k'= 1  /= 1  /'= 1 


(36) 


where  K  is  the  intensity  operator  of  the  volumetric  fluctuations.  Employing  the  weak 
fluctuation  approximation  we  approximate  K  by  its  leading  term 


□ 


>Tf)4Ai 


(37) 


On  substituting  this  in  (36)  we  get, 


E,  ®E'). =(eA .  ®(e;>,.  ■ +z  a  (g.*>,  ®(g;>,.(t<-  ®t;)(e,  ®h 


r-=i 


Next  we  average  (38)  w.r.t.  the  surface  fluctuations  to  get 


j  j 


N 


e , ® e j )  =((e,. )  ®(e;.)  )  +£  pM Gjk)  ®(g;,.)  )  (Tv®m(Et®E; 


k= 1 


where  we  have  used  the  following  approximation 

G >.  ® (G-,  >,  (T,  ®  t;  )  (E,  ®  e;  )„)  □  ((Gj,  \  ® (g;  })  {%  ®  t;  }(e,  ®  E;  \ 


(38) 


(39) 


(40) 


We  call  this  the  ‘weak  surface  correlation’  approximation,  which  we  will  see  later  to  be  an 
important  condition  embedded  in  the  RT  approach  to  our  problem. 

As  it  stands  this  equation  is  very  difficult  to  solve  either  analytically  or  numerically. 
Besides,  one  important  goal  for  us  is  to  investigate  the  conditions  needed  for  employing  the 
radiative  transfer  approach  for  our  problem.  With  these  in  mind  we  introduce  Wigner 
transforms.  Note  that  (39)  is  an  equation  for  the  coherence  function  which  is  a  ‘space- 
space’  quantity.  On  the  other  hand  the  RT  equation,  as  we  saw  earlier,  is  an  equation  for 
the  specific  intensity  which  is  a  ‘phase-space’  quantity.  Wigner  transforms  serve  as  a 
bridge  to  link  these  two  quantities  [43,9,20,30]. 
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We  introduce  Wigner  transforms  of  waves  and  Green's  functions  as 


f  y  +  r' 


,  k  =  J  (E.  (r  ))y  0  (El  ( r'))yi  exp  [-ik  ■  (r  -  r')]  d(r  -  r’) 


(41a) 


f  y  +  r' 


’k  )  =  1  (<E™  (r)>,  ®{K  (r'))X  exP \~ik  ‘  (r  _ r')] d(j  - r') 


(41b) 


^,k  [  =  J‘ d{r  “r')j d (ri  -r/) e“'k(r”r)  e'Kr'  r;)  ((Gm„(r,ri))v  ®  (G^r'.r/)^ 


(42) 


In  terms  of  these  transforms  (39)  becomes 


A 

K  (r>k)  =  K  (r’k)  +  TT^rZ  A,  j  dr'\ da\ dWnn  (r,k  I  r',a)T>',«  |  y”  ,  0)$n(r” ,  0)  (43) 

yin)  „=i  n 


with 


T,(R1,a|R2,P)  =  J  dYx  J  dv2  exp  {-/a  •  rt  +  z'P  •  r2 }  /t  (Rj  +  \ ,  R2  +  ^)®T,(R1-iR2-f) 


(44) 


where  T  is  the  element  transition  operator  in  «-th  layer. 


The  fact  that  our  problem  has  translational  invariance  in  azimuth  implies  the  following: 

=  <45a) 


(r,  k  I  r 1 )  =  sm „  ( 2,  k  |  z\  1;  v1  -  r[ ) 
Using  these  relations  in  (43)  we  have 


(45b) 


13 

Approved  for  public  release;  distribution  unlimited. 


K(z&)  =  K  (z’ k)  +  T^rrZ  A, \y  dz'j da J d$Smn  (z, k  |  z',  a;  O^F, (a,  (5)6 n (z, (5) 

(2/T)  „= i 


(46) 


where 

dmn  (*,k|  z\ a; 0)  =  | Smn  (z, k  |  z', a; r±  - r[)d(r±  - r' ) 


(47) 


#„(a,/?)  =  f(a,/?)0f*(a,/?)  (48) 

f  is  the  element  scattering  matrix  in  the  n-th  layer.  Since  the  medium  is  assumed  to  be 
sparse  inter-particle  scattering  takes  place  in  the  far-held  zone  of  each  other.  It  is  based  on 
this  fact  that  we  have  transitioned  from  T  ,  to  6n .  Also  we  have  employed  the  on-shell 

approximation  to  Tn . 

To  proceed  further  we  need  to  evaluate  Smn .  Note  that  we  need  to  relate  this  system  with 
that  of  RT,  which  involves  the  boundary  conditions  at  the  interfaces.  Therefore  we  need  to 
identify  the  coherence  functions  corresponding  to  up-  and  down-going  wave  functions.  To 
facilitate  this  we  decompose  (Gm;i)  into  its  components, 

(G  )  =8  G°  +  G"“  +G'"/  +Gdu  +Gdd  (49) 

where  the  first  term  is  the  singular  part  of  the  Green's  function.  The  superscripts  u  and  d 
indicate  up-  and  down-going  elements  of  the  waves.  The  other  components  are  due  to 
reflections  from  boundaries.  These  are  formally  constructed  using  the  concept  of  surface 
scattering  operators  as  follows  [33], 


J (7k±J dk'±  {Sf  (k1,k;)}/'1  exp{/k±  -r  +  m^(k±)z}exp{-ik'i  -r' -ibqv„(k'±)z'} 

where  is  the  surface  scattering  operator.  The  superscripts  a  and  b  on  S  are  used  to 
indicate  whether  the  waves  are  up-going  or  down-going.  In  the  exponents,  a,  b= 1  if  the 
waves  are  up-going.  We  let  a,  b  =  -1  if  the  waves  are  down-going.  The  z-component  of  the 
mean  propagation  constants  in  the  n- th  layer  is  denoted  as  qn .  We  recall  that  &mn  is  the 

Wigner  transform  of  U Gmn  'j  0  (G*nm  ^  \  .  The  superscripts  //,  v  stand  for  polarization, 
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either  h  or  v.  When  we  use  (34)  to  perform  the  Wigner  transform  we  ignore  all  cross  terms. 
In  other  words,  we  make  the  following  approximation. 


S  =8  +3*'  (51) 

mn  mn  m  mn  mn  mn  mn 

where  §„„  is  the  Wigner  transform  of  0(G “**}J  . 


With  the  introduction  of  this  representation  for  3mn  in  (43)  we  can  trace  up-  and  down¬ 
going  waves  to  obtain  the  following  equations  for  the  coherence  function: 


C'(z,k)  =  C(z,k) 


+  - 


(2k) 
4  n 


-p,n[z  dz'\ d(l] m:  (ZM  z', «; 0)T““ (a,  f3)&^  (z, P) 


(2x)4 


£p,jy  dz'\ da\ d^r  (z,k  I  z',a;0)Yamb(a,j3)£nb(z',j3) 

n= 1 


(52a) 


&m  (Z’k)  =  (Z’k) 


(2  n)‘ 
4  n 


-  pm  dz'\  da\  d&£  (z,  k  |  z',  a;  Ojrj  (a,  P)Sam  (z\  j3 ) 


(2tt)4 


£  Pn  J  '  dz'\ da\  d&Z  (z,  k  I  z',  a;  0)Tf  (a,  0)$* (z',  fi) 


(52b) 


Note  that  summation  over  a,b  =  | u ,  d  |  is  implied  in  the  above  equations.  The  first  tenn  in 


these  equations,  &sa ,  represents  the  contribution  due  exclusively  to  surface  scattering,  and 
has  the  following  form: 


n] 

J  rexP 

ia(q£-qvm  )z_ 

r“’{sa*}v 

(  m  J  (  m  J 


(k±>k±,. 


E  ,E 

pi  vi 


(53) 


where  E"  is  the  amplitude  of  the  up-going  wave  in  the  /«-th  layer  after  volumetric 
averaging  is  perfonned.  This  means  that  it  is  a  random  function  of  surface  fluctuations. 
When  we  substitute  (53)  and  the  expressions  for  Smn  in  (52)  we  find  that 


K(Z^)\m=27lS\k,  -  ~  a  (kJ_)  +  d'm  (kl )]}  eXP  \}a  (  Qm  ~  <£  )■ Z]  {K  (Zv kl)}^  (54) 


15 

Approved  for  public  release;  distribution  unlimited. 


On  substituting  this  in  (52)  and  differentiating  w.r.t.  z  we  obtain  the  following  transport  equations: 


—  ~  1  Vqu  (kl )  -  4*  (K  )]|  (Z>  M  = 

=  4  Pm  |  das;  (a± )  ®  5r  K ) '  (55a) 

•  jk± ’\[Vm (ki)  +  (ki)] ; «±4 ' a  fc' (a^)  +  & (ar )] ( z » «±) 


r/z 


i  [qu  (kj_)  -  ql  (kj]  j  <  (Z,  k±)  = 

=  4pm\da.S;(a.±)®S::  (u±)  ■ 

■  |k± » - 1  (k±)  +  ^  (k± )] ;  o± » |  a  [^-  (o± )  +  ^-  («±  )]|  K'V'  (z’ 1 °  J. ) 


(55b) 


where  summation  over  a  is  implied.  When  the  superscript  a  corresponds  to  a  the  value  of  a 
in  the  argument  of  &m  take  the  value  + 1 ;  on  the  other  hand  when  the  superscript  a 

corresponds  to  d  the  value  of  a  in  the  argument  of  JF  take  the  value  -1.  Since  all  quantities 

in  (55)  correspond  to  the  same  layer  m  we  have  dropped  the  subscript  in  in  #  and  6  to 
avoid  cumbersome  notations.  To  obtain  appropriate  boundary  conditions  we  have  to  go 
back  to  the  integral  equation  representations  for  Su  and  S*v ,  examine  their  behavior  at  the 

interfaces,  and  seek  a  relation  between  them.  After  considerable  effort  we  managed  to 
arrive  at  the  following  boundary  conditions.  At  the  (w-l)-th  interface  we  have 


f„V„-i.kJ  =  J  (ki ,  kl  >K'  (z»-.  •  kl>  *1 


(56) 


with  91  =  R  ®  R  where  R  ,  is  the  stack  reflection  matrix  (not  the  local  reflection 

m—l,m  x 

matrix)  for  a  wave  incident  from  below  on  the  (m  - 1)  -th  interface.  Similarly 
K (z„_, . k  J  =  | (k± , kl ))S!,  (z„_, , (57) 


where  91,n+1  m  is  the  tensor  product  of  stack  reflection  matrix  for  a  wave  incident  from 
above  on  the  (m  - 1)  -th  interface.  We  were  able  to  obtain  the  boundary  conditions  only 
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after  imposing  certain  approximations  such  as  the  one  given  below.  Consider  the  following 
identity: 


Qidu 
^  mm 


=  DR 


+  S 


mt 

mm 


D 


(58) 


where  Dm  =  diag  { cxp(iqlidm ),  cxp(ic/ydm )}  .  Notice  that  this  is  an  operator  relation  where  all 

elements  are  operators.  Taking  the  tensor  product  of  (58)  with  its  complex  conjugate  we 
have 

si  ®s:; = (B„  ®d:)(b„_1„  0 

{s;+sij®{s:+si}‘ 

Next  we  average  (59)  w.r.t.  surface  fluctuations  and  get 

(sdu  ®Srf“*\n(D  0D*)(R  ,  05*  ,  \ 

\  mm  mm  /  V  m  ms  \  m—\,m  m-\,m  / 

/(^  +S““  1®^  +S““  TVd  0D* ) 

\  (  m  mm  J  (_  m  mm  J  M  mm) 

where  we  have  approximated  that  the  two  tensor  products  in  the  middle  are  weakly 
correlated.  A  further  approximation  that  we  impose  is  given  as  follows: 


(D  ®D*  ) 

\  m  m  J 


(59) 


(s>  +su"  l®(s>  +  s“"  Y)n  s>  ®s>* +/S1'"  0S“"*\ 

y  m  mm  J  (_  m  mm  J  /  m  m  \  mm  mm  I 


These  are  the  kinds  of  approximations  required  to  arrive  at  our  boundary  conditions. 


(61) 
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3.4  Transition  to  Radiative  Transfer 

Next  we  have  to  transition  from  this  transport  equation  (55)  to  the  phenomenological 
radiative  transfer  equation  discussed  earlier.  To  accomplish  this  we  have  to  link  the  key 
quantities  of  waves  and  radiative  transfer,  viz.,  coherence  function  and  specific  intensity. 
The  relation  between  them  is  obtained  by  computing  the  Poynting  vector  using  the  two 
concepts.  One  of  the  fundamental  results  of  electromagnetic  wave  theory  is  the  Poynting 
vector  given  as 

(S(r)>,=iRe{(ExH-)} 

For  time-harmonic  tranverse  electromagnetic  waves  this  becomes 
<S(r)>,  =  *{(|£-„(r)|2)  +  )} 

From  Wigner  transform  relations  we  have 

(£„  (r)£„  (r)>  =  ^  J S(„  (r,  kjrfk,  (62) 

The  average  Poynting  vector  is  also  related  to  the  specific  intensity  as 

(S(r))  =  |M(r,k>/Q  (63) 


The  above  two  relations  suggest  that  the  following  definition  for  the  specific  intensities 


I„v(z>s) 


1  k'2 

2 77  (2/r)2 


cos^v(z,k±) 


(64) 


Now  we  can  transition  to  the  phenomenological  RT  equations.  Using  the  relation  between 
6  and  /  we  change  the  integration  variable  to  solid  angle  and  arrive  at  the  following 
equations, 


z)  d 

COS# - 1-/7.. 

dz  h} 


■  =  (n, n')/“(z, »')+/;"  (n,n')/'(z,i')pn'  (65a) 


a  d  , 

-cos# - hr?.. 

dz  J 


' (z,s)  =  [p.,D!)rj{z,s')  + Pyd  (Q,Q')/"'(z,,?')}#Q' 


(65b) 


where  77  is  the  extinction  matrix  and  P  is  the  phase  matrix.  Implicit  summation  over 
subscript  j  is  assumed  in  (65).  To  facilitate  comparison  with  the  results  of  Ulaby  et  al  [38] 
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and  Lam  and  Ishimaru  [14]  we  have  used  a  modified  version  of  Stokes  vector  [11].  Instead 
of  the  standard  form  {l,Q,U,V}  we  use{(/  +  0/2, (7-0/2, t/,F}  .  The  subscript  of  I  in 

(65)  denotes  the  element  number  of  our  modified  Stokes  vector.  Although  the  structure  of 
this  equation  is  identical  to  that  of  the  RT  (equation  (2)),  the  elements  of  the  phase  matrix 
and  the  extinction  matrices  are  not  the  same.  The  primary  reason  is  because  of  the 
differences  in  the  real  part  of  the  mean  propagation  constants  of  horizontally  and  vertically 
polarized  waves.  On  assuming  that  q'h  =  q'v  =  k'mz  we  obtain  the  following  expressions  for 
the  extinction  and  phase  matrices: 


>1  =  —  cos  0  diag  { 2q",  2 q"h ,  q"  +  q"h ,  q"  +  q"h } 

(66) 

p;‘(n,n')  =  g;i{ 

k ±,ak  cos  0\  k/ , bk  cos  O'] 

(67) 

where 

ii 

Fo 

%=(i/.»r) 

%=<R  e(/„/.;)) 

^4  =-<Im(/„/,;)> 

3,=(i  fj) 

rN _ 

-s; 

ii 

CN 

=  2  {Re(/„ /„;,)> 

%=2(r  <=(/,, ,/;,)) 

^3 = ^  (Re(/w/*;  >) 

+  {Re(/v„/;v)) 

«,  =  2(h(U;)) 

^=(lm  (/„/!)) 
+(Im<XX)) 

-{Re(/„/;,)) 

(68) 

We  have  suppressed  the  arguments  for  brevity.  For  instance, 
{k1,akcos<9;k'L,Mcos6''} 

=  ^Re|/W  (it  L,ak  cos  0;k'L,bk  cos  0'}  f*h  (k1,«£cos6>;k'L,bkcos6,')j^ 

/  '  s  are  the  elements  of  the  scattering  matrix  f  of  particles  defined  as  follows : 


E  =  —  f  E. 

.v  r  i 


In  the  {/7,  v}  basis  f  is  given  as 
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Note  that  these  transport  equations  (65)  are  identical  to  those  of  classical  RT  equations  (2) 
that  we  described  in  Section  2.  Thanks  to  our  statistical  wave  approach  we  now  have 
explicit  expressions  for  the  extinction  matrix  and  phase  matrix  in  terms  of  the  statistical 
parameters  of  the  problem.  Let  us  now  next  turn  our  attention  to  the  boundary  conditions 

(BC).  In  our  wave  approach  we  obtained  BCs  in  terms  of  ‘stack’  reflection  matrix  R  , 
whereas  in  the  RT  approach  the  BCs  are  given  in  terms  of  the  local  interface  reflection 
matrices.  We  can  readily  reconcile  with  this  apparent  difference.  Note  that  the  BC  in  the 
wave  approach  forms  a  closed  system  whereas  in  the  RT  approach  it  is  open  (linked  to 
adjacent  layer  intensities).  Let  us  take  a  look  at  the  BC  at  the  (m- 1  )-th  interface.  R  m  can 

be  expressed  in  terms  of  Rm_2  m_,  as  follows, 


R  =  R  +  T 

\,m  tm— l,m  1 


-1-1 


i  —  r  n  d 
1 


R  D  T 

m— 2,m-r  m- 1  \,m 


(70) 


This  is  the  relation  between  the  stack  reflection  coefficients  of  adjacent  interfaces.  The  R 
and  T  are  local  (single  interface)  reflection  and  transmission  matrices  at  the  (m  - 1)  -th 
interface.  On  operating  E"  with  (70)  we  get 


Ed  =  R  .  E"  +T  ,E 

m  m-\,m  m  m.m- 1 


d 

m- 1 


(71) 


Notice  that  this  boundary  condition  now  involves  only  local  interface  Fresnel  coefficients. 
Take  the  tensor  product  of  (71)  with  its  complex  conjugate  and  average  w.r.t.  surface 
fluctuations.  Employing  the  Wigner  transform  operator  on  this,  we  obtain  a  boundary 
condition  at  the  (m- 1  )-th  interface  similar  to  that  of  the  RT  system.  However,  the  reflection 
and  transmission  matrices  used  in  the  RT  system  correspond  to  unperturbed  medium  as 
opposed  to  the  average  medium  as  in  the  case  of  the  wave  approach. 


Similarly  we  write  R  ,  in  terms  of  R  .  ,  and  hence  obtain  the  BC  at  the  m- th 

J  m  +  \,m  m+2,m+\ 

interface  as 


E"  =  R  E  +T  E" 

m  m+\,m  m  m,m+ 1  m+ 1 


(72) 


Take  the  tensor  product  of  (72)  with  its  complex  conjugate  and  average  w.r.t.  surface 
fluctuations.  Employing  the  Wigner  transform  operator  on  this,  we  obtain  boundary 
condition  at  the  m- th  interface  identical  to  (3)  (after  making  the  approximation  as  before). 
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4.  Results  and  Discussion 


The  main  goal  of  this  report  is  to  critically  examine  the  radiative  transfer  (RT)  approach  to 
remote  sensing  of  layered  random  media  with  rough  interfaces.  Such  problems  are  often 
encountered  in  applications.  Several  assumptions  are  embedded  in  the  radiative  transfer 
approach.  There  are  numerous  works  in  the  literature  on  the  study  of  radiative  transfer 
theory  and  the  underlying  assumptions.  However,  all  have  dealt  with  unbounded  geometry 
or  bounded  geometry  with  nonscattering  boundaries.  Our  interest  in  this  report  is  on  the 
problem  of  layered  random  media  with  irregular  scattering  boundaries.  There  does  not  exist 
any  critical  study  of  the  RT  approach  to  this  kind  of  geometry.  Our  study  has  shown  that 
there  are  additional  conditions  embedded  in  the  RT  approach  to  this  problem  than  for  the 
problem  with  unbounded  or  nonscattering  geometries.  These  facts  are  not  well  known  to 
the  users  of  the  RT  approach.  One  purpose  of  this  report  is  to  inform  the  remote  sensing 
community  about  these  additional  conditions  so  that  they  have  a  good  idea  on  when  the  RT 
approach  is  an  acceptable  model  for  the  application  at  hand. 

To  enable  this  study  we  developed  a  statistical  wave  theory  for  the  combined  problem  of 
layered  random  media  with  rough  interfaces.  Such  a  foundation  is  essential  for  this  study 
because  there  is  no  suitable  wave  theory  in  the  literature  that  is  suitable  for  multiple 
scattering  for  our  problem.  From  this  study  we  find  that  the  coherent  waves  for  our 
problem  behave  like  waves  in  a  layered  homogeneous  media  with  planar  interfaces.  The 
propagation  constants  of  this  layered  media  are  primarily  determined  by  the  statistical 
properties  of  the  local  medium.  It  is  weakly  dependent  on  interface  roughness.  Its 
dependence  on  medium  fluctuations  of  other  layers  is  of  higher  order.  In  contrast,  the 
effective  reflection  coefficients  are  influenced  by  not  only  statistical  properties  of  the  local 
interface  but  also  on  the  statistical  properties  of  the  adjacent  media.  Its  dependence  on  the 
surface  roughness  of  other  interfaces  and  media  are  of  higher  order.  We  notice  that  diffuse 
scattering  is  fairly  complicated  because  of  volume-surface  interactions.  In  recent  years 
there  has  been  a  proposition  that  the  combined  problem  of  volumetric  scattering  and 
surface  scattering  be  split  into  two  parts:  one  due  to  volumetric  scattering  and  the  other 
surface  scattering.  The  surface  scattering  part  is  due  to  the  layered  structure  with  rough 
interfaces  and  homogeneus  media  with  effective  permittivities.  The  volumetric  scattering  is 
from  the  layered  random  medium  with  unperturbed  at  interfaces.  It  is  clear  from  our  results 
that  this  kind  of  splitting  is  not,  in  general,  possible  unless  the  fluctuations  of  the  problem 
are  very  weak. 

The  results  obtained  in  this  report  apply  to  two  situations.  Below  we  consider  them 
separately: 

(a)  Time- Varying  Problem:  This  means  that  the  parameters  of  the  problem  vary  with  time. 
However,  an  important  assumption  must  be  made  to  simplify  the  analysis.  The  time 
constant  associated  with  the  parameter  fluctuations  of  the  problem  should  be  much  larger 
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than  the  time  constant  associated  with  signal  of  the  incident  wave.  The  time  average  is 
taken  over  a  period  much  larger  than  the  time  constant  associated  with  the  medium 
fluctuations  (that  includes  volumetric  and  surface  fluctuations).  Now  we  impose  the 
ergodicity  hypothesis  and  equate  time-averages  to  ensemble-averages.  One  such 
application  where  this  situation  occurs  is  in  electromagnetic  wave  propagation  and 
scattering  in  atmosphere. 

(b)  Time  Independent  Problem:  Here  all  the  parameters  are  independent  of  time.  However 
they  undergo  spatial  variations.  The  statistical  aspects  of  the  problem  enter  through  spatial 
fluctuations.  Therefore  all  averages  are  ensemble  averages.  Obviously  then  the  question  of 
ergodicity  does  not  arise  here. 

The  transport  equations  as  derived  in  (55)  along  with  the  boundary  conditions  (56)  and  (57) 
are  important  results  of  the  report.  This  system  describes  the  behaviour  of  up  and 
downward  travelling  coherence  functions  in  each  layer  of  the  problem.  Several  important 
physical  quantities  can  be  directly  calculated  using  these  coherence  functions. 

(a)  Poynting-Stokes  Tensor:  The  Poynting-Stokes  tensor  [21,22],  a  key  quantity  in  radiative 
transfer,  is  related  to  our  coherence  function  as  follows: 

P(z)  =  J  k+  x S"(z,k±)t/k±  +  \  k~  x  Sd(z,ls.Jdk±  (73) 

[n+  cr  J 


where  ij{z)  is  the  intrinsic  impedance  of  the  medium  where  z  is  located,  k 1  and  k  are  the 


unit  propagation  vectors  associated  with  &u  and  Sd  respectively.  Q  and  Q  indicate  that 
the  domain  of  integration  is  the  upper  and  lower  hemisphere  respectively. 

(b)  Specific  Intensity:  We  showed  the  relationship  between  coherence  function  and  the 
specific  intensity  governed  by  the  phenomenological  radiative  transfer  equation  (65). 
Details  are  given  below: 


x~ 

&yv 

/ 
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—  i  ^  cos  0 

^  (2 4 
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V  _ 
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(c)  Average  Poynting  Vector:  This  quantity  which  is  vital  for  radiation  budget  computation 
[36,18]  is  derived  from  the  coherence  function  as  follows: 

{s(z)>=*te7iZ  |  jdk1  +  j  rs;(z.k1)aj  (75) 

^  [n+  cr  J 


(d)  Bistatic  Scattering  Coefficient:  In  our  problem  the  system  is  excited  by  a  plane  wave 
incident  on  the  zero-th  interface  from  above.  The  bistatic  scattering  coefficient,  cr'  in 
Region  0  (in  the  direction  k±  )  is  given  as 
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<r'(ki.ku) -afp^rj  (^(ki.kD^-fz.kDAl  (76) 

where  9C01  =  Xm  ®  X*n  .  Suppose  the  incident  plane  wave  is  horizontally  polarized.  Then  the 
lih  component  of  LHS  represents  cr'hh ;  The  vv  component  represents  <jrvh .  Similarly  if  the 
incident  wave  is  vertically  polarized,  the  hh  component  of  RHS  represents  cr'hv ;  the  vv 


component  represents  arm .  We  have  used  the  superscript  r  to  indicate  that  this  quantity  is  a 


reflection  type  scattering  coefficient  since  the  source  and  observation  points  are  located  in 
the  Same  Region  0.  This  kind  of  scattering  coefficient  is  extensively  used  in  remote  sensing 
applications  [38,8,1,37].  Next  suppose  that  the  observation  point  is  in  Region  N  +  \  while 
the  source  is  still  in  Region  1 .  Now  the  bistatic  scattering  coefficient  is  given  as 


cr'fk 


i»M  = 


cos# 


cos0i  (2*y 


■|(w  k^kD^z^ki)*! 


(77) 


We  call  this  the  transmission  scattering  coefficient  and  indicate  it  with  the  superscript  t. 
This  quantity  is  quite  useful  in  a  variety  of  imaging  applications. 

(e)  Passive  Remote  Sensing:  The  key  quantity  of  interest  in  passive  remote  sensing  is 
emissivity  e  [8,38,1]  given  as 

eP=l~rpc-rPi  P  =  {h,v]  (78) 


where  r  ,  and  r  .  are  the  coherent  and  incoherent  reflectivities  defined  as  shown  below. 

pc  pi 


V=Z|(4(ki,k1,))[ 

<1 

'V  =Zj</(kL’k^>/k 


(79a) 

(79b) 


R00  is  the  stack  reflection  coefficient.  <jr  is  the  reflection  scattering  coefficient  defined 
earlier. 


Indeed  there  are  other  useful  physical  quantities  [8,19,13]  that  can  be  obtained  from  the 
coherence  function  derived  in  (55).  We  have  just  given  some  examples  relevant  to  remote 
sensing  problems.  In  order  to  understand  the  foundations  of  the  radiative  transfer  approach 
we  made  the  transition  from  the  governing  equation  for  coherence  function  to  the 
phenomenological  radiative  transfer  equations.  We  made  the  transition  from  statistical 
wave  theory  to  radiative  transfer  theory  by  employing  the  Wigner  transfonn  and  computing 
the  average  Poynting  vector  by  using  wave  theoretical  methods  and  transport  theoretical 
concepts. 

Having  made  the  transition  from  statistical  wave  theory  to  radiative  transfer  theory  it  is 
now  instructive  to  itemize  the  assumptions.  The  three  basic  conditions  required  are: 

I.  Quasi-stationary  field  approximation 


23 

Approved  for  public  release;  distribution  unlimited. 


II.  Sparse  distribution  of  scatterers 

III.  The  number  of  particles  in  each  layer  is  large 

These  are  the  well-known  conditions  necessary  for  the  unbounded  random  medium 
problem.  However,  if  the  medium  is  bounded  we  need  to  impose  additional  conditions.  We 
found  that  the  extinction  coefficients  calculated  in  the  wave  approach  and  the  RT  approach 
are  different  and  only  after  further  approximations  can  they  be  made  to  agree  with  each 
other.  The  additional  condition  required  for  our  bounded  random  medium  problem  is: 

IV.  Layer  thickness  must  be  of  the  same  order  or  greater  than  the  mean  free  path  of  the 
layer. 

When  the  interfaces  are  randomly  rough  we  further  need  the  following  conditions. 

V.  Weak  surface  correlation  approximation. 

VI.  All  fluctuations  of  the  problem  are  statistically  independent  of  each  other. 

The  above  are  the  essential  conditions  associated  with  the  RT  approach.  There  are  several 
secondary  assumptions  that  we  have  employed  for  simplifying  the  analysis  and  discussion. 
They  are  given  as  follows: 

1.  The  location  and  orientation  are  the  only  randomly  varying  aspects  of  the  media. 

2.  All  particles  in  a  particular  layer  are  of  the  same  type  and  are  uniformly  distributed. 
However  they  can  be  of  different  type  in  different  layers. 

3.  The  problem  is  translationally  invariant  and  isotropic  in  the  azimuth. 

4.  The  problem  is  assumed  to  be  time-independent. 

5.  If  the  problem  happens  to  be  time-varying  the  following  additional  conditions  need  to  be 
imposed:  (a)  The  time  constant  of  statistical  parameters  is  much  longer  than  that  of  the 
signal,  (b)  The  problem  is  ergodic. 

6.  The  interfaces  of  the  random  medium  are  parallel  planes  on  the  average. 

Let  us  make  some  remarks  on  the  secondary  assumptions. 

Assumption  1 :  One  can  add  to  it  other  variations  such  as  size,  material  properties,  etc.,  at 
the  price  of  more  complicated  analysis. 

Assumption  2:  All  particles  need  not  be  of  the  same  type.  We  can  have  multi-species  in 
each  layer.  The  distributions  can  be  arbitrary.  It  is  not  necessary  for  them  to  be  uniform. 
However,  such  details  will  add  to  the  complexity  of  the  results. 

Assumption  3:  While  translational  invariance  is  an  important  assumption  for  our  analysis, 
the  statistical  fluctuations  need  not  be  isotropic  in  azimuth.  Layered  random  media 
without  translational  invariance  is  a  very  complex  problem  not  considered  in  this 
report. 

We  can  indeed  take  into  account  anisotropy  in  statistical  fluctuations.  However,  we  will 
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have  a  tensor  form  of  effective  permittivity  and  associated  mode  patterns. 

Assumption  4\  Time-independence  is  another  assumption  intended  to  simplify  the  analysis. 
Assumption  5:  The  results  obtained  in  the  report  apply  to  time-varying  case  provided 
additional  assumptions  are  imposed. 

Assumption  6:  This  is  an  important  assumption  necessary  for  our  analysis.  However  this 
assumption  is  not  essential  for  the  main  character  of  our  results. 

We  would  like  to  make  a  couple  of  more  remarks  before  ending  the  discussion  of  our 
assumptions:  (a)  In  RT  theory  the  medium  is  assumed  to  be  sparse  and  hence  the 
“refraction  effects”  of  the  fluctuations  are  ignored.  Thus  in  the  boundary  conditions  we 
should  use  the  background  medium  parameters  rather  than  the  effective  medium 
parameters  as  derived  in  our  statistical  wave  theory,  (b)  To  arrive  at  (65)  we  have  ignored 
the  contribution  of  evanescent  modes. 

Recently  Mishchenko  et  al.  [24]  (hereafter  referred  to  as  MTL  for  brevity)  derived  the 
vector  radiative  transfer  equation  (VRTE)  for  a  bounded  discrete  random  medium  using  a 
rigorous  microphysical  approach.  This  enabled  them  to  identify  the  following  assumptions 
embedded  in  the  VRTE. 

1 .  Scattering  medium  is  illuminated  by  a  plane  wave. 

2.  Each  particle  is  located  in  the  far  field  zone  of  all  other  particles  and  the  observation 
point  is  also  located  in  the  far  field  zones  of  all  the  particles  forming  the  scattering 
medium. 

3.  Neglect  all  scattering  paths  going  through  a  particle  two  or  more  times  (Twersky 
approximation). 

4.  Assume  that  the  scattering  system  is  ergodic  and  averaging  over  time  can  be 
replaced  by  averaging  over  particle  positions  and  states. 

5.  Assume  that  (i)  the  position  and  state  of  each  particle  are  statistically  independent 
of  each  other  and  those  of  all  other  particles  and  (ii)  spatial  distribution  of  the 
particles  throughout  the  medium  is  random  and  statistically  uniform. 

6  Assume  that  the  scattering  medium  is  convex. 

7.  Assume  that  the  number  of  particles  N  forming  the  scattering  medium  is  very  large. 

8.  Ignore  all  the  diagrams  with  crossing  connections  in  the  diagrammatic  expansion  of 
the  coherency  dyadic. 

First,  notice  that  the  MTL  dealt  with  the  problem  having  nonscattering  boundary. 

For  this  problem,  only  the  first  three  of  our  essential  conditions,  and  the  first  five  of  our 
secondary  conditions  should  apply.  The  list  of  assumptions  given  by  MTL  consists  of  a 
combination  of  both  essential  and  secondary  conditions.  Furthermore,  the  conditions  that 
MTL  obtained  are  not  identical  with  the  ones  that  we  derived  even  when  restricting 
attention  to  the  random  medium  problem  with  nonscattering  boundary.  This  is  because  the 
methods  we  employed  are  different.  Here  below  we  address  the  differences  and  provide 
explanation  for  the  discrepancies. 
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MTL  2:  MTL  have  made  deliberate  use  of  the  far  field  approximation.  We  have  not 
explicitly  used  this  approximation.  Instead  we  have  used  the  more  general  quasi¬ 
uniform  field  approximation.  To  illustrate  this  consider  electric  fields  at  two  points 
r,  and  r2 : 

E(ii)  =  A(riykr'  E(r2)  =  A(r2ykr2  (80) 

The  coherence  function  of  these  fields  is  given  as 

(E(r1)®E*(r2))  =  (A(r1)®A*(r2))e'k(r-r2)  (81) 

In  quasi-stationary  field  approximation  we  assume  that 

(A(r, )  ®  A*(r2 ))  =  <?((r,  +  r,)  /  2)  (82) 

where  dP  :=  ^A®  A*^ .  In  contrast,  MTL  employed  the  far-field  approximation  under 

which  A  is  independent  of  position  vectors.  For  our  problem  of  layered  random  media 
explicit  implementation  of  far  field  approximation  is  unrealistic. 

MTL  3:  We  have  not  explicitly  used  the  Twersky  approximation.  Instead  we  made  use  of 
the  fact  that  the  number  of  particles  is  large  and  hence  approximated  conditional 
averages  by  unconditional  averages.  Incidentally,  we  arrive  at  the  same  result  as  that  of 
MTL  for  the  unbounded  random-  medium  problem. 


MTL  6:  In  our  problem  we  have  distinct  scattering  boundaries  and  the  character  of  the 
waves  exiting  or  entering  them  are  explicitly  contained  in  the  boundary  conditions. 
Hence  convexity  of  the  scattering  medium  is  not  a  necessary  condition  for  us. 

MTL  8:  In  our  approach  we  used  the  weak  fluctuation  criterion  and  retained  only  the 
leading  term  of  the  intensity  operator.  Ladder  term  is  the  leading  tenn.  Crossed  terms 
are  higher  order  terms  which  have  contribution  only  in  the  backscattering  direction. 
Hence  such  terms  do  not  appear  in  our  results. 

Let  us  reiterate  that  MTL  [24]  with  whom  we  made  detailed  comparison  have  considered  a 
bounded  random  medium  with  nonscattering  boundary.  Ours  is  a  significantly  more 
complicated  problem  because  of  the  layered  structure  and  scattering  interfaces  with 
irregular  geometry.  This  means  that  our  problem  contains  volume  surface  interactions 
because  of  multiple  scattering.  Nobody  has  critically  examined  the  RT  approach  to  this 
problem.  By  following  a  systematic  and  rigorous  approach  we  found  that  the  RT  approach 
to  our  problem  requires  additional  conditions  beyond  those  required  for  the  problem  with 
nonscattering  boundaries.  It  is  expected  that  some  of  our  conditions  are  equivalent  to  those 
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of  MTL.  The  key  results  of  the  report  lie  in  the  extra  conditions.  These  additional 
conditions  that  are  special  to  our  geometry  are:  (a)  The  layer  thickness  must  be  of  the  same 
order  or  larger  than  the  mean  free  path  of  the  medium  (b)  Weak  surface  correlation 
approximation,  (c)  All  volumetric  fluctuations  and  surface  fluctuations  are  independent  of 
each  other.  These  are  the  essential  extra  conditions.  Following  are  some  secondary  extra 
conditions  that  we  have  imposed  to  arrive  at  our  results.  The  surface  fluctuations  are:  (a) 
translationally  invariant  and  isotropic  in  azimuth,  (b)  statistically  homogeneous,  and  (c) 
single-valued.  These  secondary  assumptions  were  employed  in  order  to  simplify  the 
analysis  and  discussion. 

The  existence  of  these  extra  conditions  means  that  RT  approach,  as  it  is  popularly 
conceived,  may  be  applied  to  only  a  limited  type  of  layered  random  media.  In  other  words, 
the  accuracy  which  can  be  obtained  using  RT  theory  depends  on  the  geometry  of  the 
problem.  For  certain  geometries  such  as  the  one  discussed  in  this  report  the  results  obtained 
using  RT  approach  can  be  grossly  inaccurate.  We  have  obtained  these  results  using  a 
systematic  analysis  of  the  macroscopic  Maxwell's  equations.  Nobody  else  has  observed  that 
additional  conditions  are  involved  in  the  RT  approach  to  layer  geometry  than  for  the 
unbounded  geometry  of  with  nonscattering  boundaries.  In  view  of  these  remarks,  the 
results  of  this  report  are  important  to  the  remote  sensing  community.  Until  now,  RT  theory 
has  been  taken  as  a  fundamental  law  and  applied  to  variety  of  random  media  with  complex 
geometries.  This  report  shows  that  under  what  conditions  RT  theory  may  be  used  to  layered 
random  media. 
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5.  Conclusion 


To  summarize,  we  have  enquired  into  the  assumptions  involved  in  adopting  the  radiative 
transfer  approach  to  scattering  from  layered  random  media  with  rough  interfaces.  To 
facilitate  this  enquiry  we  adopted  a  wave  approach  to  this  problem  and  derived  the 
governing  equations  for  the  first  and  second  moments  of  the  wave  fields.  We  employed 
Wigner  transforms  and  transitioned  to  the  system  corresponding  to  that  of  radiative  transfer 
approach.  In  this  process  we  found  that  there  are  more  conditions  implicitly  involved  in  the 
RT  approach  to  this  problem  than  it  is  widely  believed  to  be  sufficient.  With  the  recent 
development  of  fast  and  efficient  algorithms  for  scattering  computations  and  the  enormous 
increase  in  computer  resources  it  is  now  feasible  to  take  an  entirely  numerical  approach  to 
this  problem  without  imposing  any  approximations.  In  spite  of  such  developments,  to  keep 
the  size  of  the  problem  manageable  only  special  cases  have  been  studied  thus  far 
[10,31,29,33].  Hence  it  is  very  much  of  relevance,  interest  and  convenience  to  apply  the  RT 
approach  to  these  problems.  However,  one  should  keep  in  mind  the  assumptions  involved 
in  such  an  approach.  Otherwise  interpretations  of  results  based  on  RT  theory  can  be 
misleading. 
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